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The plaquette phase of the square lattice quantum dimer model is studied using a continuous-time 
reptation quantum Monte Carlo method for lattices of sizes up to 48 x 48 sites. We determine the 
location of the phase transition between the columnar and plaquette phases to occur at V c /J = 
0.60 ±0.05 which is significantly larger than inferred from previous exact diagonalization studies on 
smaller lattices. Offdiagonal correlation functions are obtained. They exhibit long-range order in 
the plaquette phase but not at the Rokhsar-Kivelson point. We also observe significant finite-size 
corrections to scaling for the transition between the plaquette phase and the critical resonating 
valence bond liquid. This study demonstrates the importance of understanding finite-size effects 
when considering critical properties of the square lattice quantum dimer model. 
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I. INTRODUCTION 

The Quantum dimer models(QDMs) are interesting 
as they constitute simple examples of effective quantum 
lattice models with restricted Hilbert spaces. Quantum 
lattice models are frequently encountered in condensed- 
matter physics. They are defined by an Hamiltonian act- 
ing on a Hilbert space which is a direct product of Hilbert 
spaces for each site. 

In many cases the Hamiltonian consists of terms of 
widely different magnitude. The effects of the largest 
terms on the low-energy physics can then be effectively 
described by restricting the allowed Hilbert space. Some- 
times this can be carried out at the level of the site 
Hilbert spaces, just restricting the allowed states on a 
single site. However in other cases the constraints cannot 
be modelled in this site-specific way and the constraint 
on the Hilbert space is more complicated. The QDMs 
are examples of the latter. 

The QDMs were originally proposed* as models for an- 
tiferromagnets where the tendency to form short-range 
spin-singlet valence bonds is strong. This formation of 
short-range valence bonds is modelled as a constraint 
which effectively reduces the Hilbert space to a nontrivial 
state space; that of dimer-coverings of the lattice^. 

While the Hilbert space of a QDM is special its Hamil- 
tonian is rather simple. It consists of a kinetic term that 
flips the orientation of parallel dimers and a potential 
energy that associates an energy cost/gain to parallel 
dimers. The QDMs can be formulated on any lattice and 
their phase diagrams are largely similar. They consist 
of a liquid and various solid phases. The liquid phase is 
known as the resonating valence bond(RVB) liquid. The 
solid phases comes in at least three varieties, one phase 
with a maximum amount of parallel dimers; the columnar 
phase, one with no parallel dimers; the staggered phase, 
and one intermediate phase characterized by having a set 
of plaquettes with parallel dimers changing orientation 
constantly; a resonating plaquette phase. 

While there are good evidences for the existence of 
a columnar to a plaquette state phase transition in the 



QDM on the hexagonal 3 and triangular- lattices, the ev- 
idence is weaker on the square lattice where it is only 
based on exact diagonalization studies of linear system 
sizes up to L = In this article we show evidence of 
the columnar to plaquette phase transition for the square 
lattice QDM and estimate its location. We find that the 
plaquette phase is realized in a much smaller region in 
parameter space than previously estimated. 

The resonating plaquette phase is most directly char- 
acterized by order in a quantity offdiagonal in the dimer 
basis that measures resonating dimers. We have therefore 
carried out simulations focusing on the possible appear- 
ance of long-range order in the offdiagonal dimer flip cor- 
relation function. While we find long-range-order in this 
quantity inside the plaquette phase, the magnitude of the 
order is significantly reduced from the value expected for 
an ideal plaquette product state with resonating dimers 
on one of the four plaquette sublattices, see Fig. 2] 

The ground state of a QDM is explicitly known at its 
Rokhsar-Kivelson(RK) point 6 .. The particular form of 
the ground state implies that any T=0 quantum corre- 
lation function of observables that are diagonal in the 
dimer basis can be obtained as an infinite temperature 
correlation function of the classical dimer model, i.e they 
are just properties of the dimer-coverings themselves. 
This mapping has been utilized to calculate ground-state 
properties at the RK-poinliiSiS. On the square lattice the 
ground state at the RK-point is critical. It was initially 
considered likely that the properties of the square lat- 
tice RK-point would extend also to the immediate vicin- 
ity of the RK-point 6 . However numerical diagonalization 
studies concluded that a significant portion of the phase 
diagram exhibits crystalline order«Afi and that there were 
no evidence for a liquid state away from the RK-point s . 
These studies were carried out on rather small lattices 
(up to L = 8) and properties of the corresponding phase 
transition at the RK-point were not addressed. In this ar- 
ticle we also attempt to address these properties by mea- 
suring the Binder ratio of the columnar order parameter 
close to the RK-point for system sizes up to L = 48. 

The quantum Monte Carlo method employed in this 
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FIG. 1: Schematic drawings of the ideal columnar state (left) 
and the ideal resonating plaquette state (right). The thick 
lines represent a high average dimer probability (1 for the 
ideal columnar state and 1/2 for the ideal plaquette state). 
The letters in the right panel show the assignment of the four 
different plaquette sublattices. 



article is a synthesis of the Continuous-time lattice Dif- 
fusion Monte Carlo method introduced in Refpii and 
the Reptation Monte Carlo method introduced in Ref>i£. 
This amalgam of methods which we term Continuous- 
Time Reptation Monte Carlo (CTRMC) is easy to 
implement and is free of population-bias and time- 
discretization errors that hamper various other forms of 
projector Monte Carlo techniques. The method is not 
restricted to QDMs and can be applied to any quantum 
lattice model free of the sign problem. 

Before we explain our results we will discuss the phase 
diagram of the square lattice QDM in greater detail. 




FIG. 2: Color online. Energy per plaquette for the lowest 
lying state in three different topological sectors characterized 
by the transition graph winding numbers (w x ,w y ) with re- 
spect to the columnar reference configuration. From top to 
bottom the energy curves have winding numbers (1, 1),(1,0) 
and (0,0). 
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II. PHASE DIAGRAM 



The square lattice QDM Hamiltonian is 
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where the summations are taken over all elementary pla- 
quettes of the lattice. We will choose units of energy such 
that the flipping energy J = 1. 

The state space of the square lattice QDM is naturally 
divided into separate topological sectors each invariant 
under the action of the Hamiltonian. Any dimer config- 
uration belongs to a topological sector characterized by 
the winding numbers of its transition graph to a refer- 
ence configuration, which we take to be the ideal colum- 
nar state shown in Fig.^left panel. The transition graph 
is obtained by overlaying the reference configuration on 
the dimer configuration in question and erasing overlap- 
ping dimers. This leaves a set of loops which might wind 
around the lattice. For V < 1 the topological sector with 
zero winding numbers have the lowest energy, see Fig. [21 
and we will restrict our simulations to this topological 
sector. 

A schematic zero temperature phase diagram of the 
QDMs is shown in Fig. El For V = 1, the RK-point, the 
ground state is the equal-amplitude superposition of all 
dimer coverings of the lattice. For the square lattice this 



FIG. 3: Generic T — phase diagram of QDMs. The phases 
are from left: columnar phase, resonant plaquette phase, RVB 
liquid, staggered phase. For QDMs on bipartite lattices: V p = 
1. 



implies that dimer-dimer correlations have no long-range 
order, but are critical and decay as a power lawi We 
will refer to this state as the RVB liquid although it is 
gapless and is believed to exist only at a single point in 
the phase diagram for QDMs on bipartite latticesiS*^. 
For non-bipartite lattices like the triangular^ and the 
Kagome lattice^ the RVB liquid has gapped excitations 
and extends over a finite region in parameter space. 

The QDMs exhibits a number of crystalline phases. 
For V — > — oo the ideal columnar state with a maximal 
number of parallel dimers will be preferred, see Fig. ^ 
This four-fold degenerate (on a square lattice) state is di- 
agonal in the dimer basis thus the kinetic term will tend 
to destroy it. These quantum fluctuations will for finite V 
lead to "disorder" within the columns, but it is expected 
that the broken rotational symmetry of the ideal colum- 
nar state still survives at least up to a critical value of V. 
Thermal effects will also tend to destroy the V — > — oo 
columnar state, these were studied in RefJ^. It is rea- 
sonable to believe that the kinetic term will eventually, 
for big enough values of V, turn the ground state into 
a state resembling the ideal resonating plaquette state 
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shown in Fig. ^ The ideal resonating plaquette state is 
also four-fold degenerate on a square lattice and can be 
written as 



IV'plaq) 



(2) 



where the product is taken over all plaquettes on one of 
the four plaquette sublattices. A similar resonant plaque- 
tte phase is known to exist for QDMs on the hexagonal 3 
and triangular— lattices. 

While the ideal columnar state is certainly the true 
ground state for V — ► — oo there is no guarantee that 
the ideal plaquette state is the ground state for V = 
as the ideal plaquette state is not an eigenstate of all 
the kinetic terms. In fact it is only an eigenstate of the 
subset of kinetic terms acting on the plaquettes with res- 
onating dimers. Nevertheless the V = point is believed 
to be situated inside the plaquette phase for both the 
hexagonal 3 -, triangular— and the square^ lattices. 

It is clear that when \V\ < 1 the ground state should 
have an appreciable amount of resonating dimers. In or- 
der to appreciate the deviations from the idealized states 
depicted in Fig.^ an d the significance of the resonating 
dimers we have plotted the average dimer densities for 
several values of V in Fig.Q] These plots and Fig. |3 were 
obtained using the Monte Carlo method which will be 
explained in the next section. In order to break transla- 
tional and rotational symmetry the plots were obtained 
on a rectangular lattice with open boundary conditions. 
For V < 1 (all panels except the bottom left in Fig. 0} 
one can see that a maximal number of flippable plaque- 
ttes is favored. This is seen as the darker squares in the 
plots. For V — —0.5 (upper left panel) the dimers on 
the central plaquette prefer to be horizontally aligned. 
We take this as an indication of the columnar phase al- 
though for the small system shown here it might as well 
be a boundary effect. The orientational preference weak- 
ens as V increases, and at V = 1, the RFC-point, the 
plaquette pattern is barely visible at all. In fact it dis- 
appears in the thermodynamic limit, and the remnant 
seen here is a finite-size effect. For V > 1 a staggered 
arrangement of dimers is preferred. This is seen to be 
true locally in the lower left panel of Fig.0] However the 
boundary conditions restrict the type of configurations 
that are available globally and the actual global pattern 
seen is a result of freezing into a specific configuration 
that depends on initial conditions and the exact pattern 
of quantum fluctuations. The V > 1 phase of the square 
lattice QDM model is also very sensitive to perturba- 
tions as shown in Refii&. We will consider V < 1 in the 
remainder of this article. 

The primary aim of this article is to find the location 
V c of the phase transition between the columnar and pla- 
quette phases. However first the quantum Monte Carlo 
method used in this article will be explained. 
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FIG. 4: Average dimer densities on a rectangular lattice with 
open boundary conditions for different values of V. The gray- 
scale marks the dimer density. Clockwise from upper left: 
V=-0.5, 0.8 ,1.0, 1.2. 



III. MONTE CARLO METHOD 



The object of CTRMC is to perform a stochastic sim- 
ulation of the imaginary-time evolution operator e~ Hr 
which in the large time limit is a projection operator onto 
the ground-state of the Hamiltonian. As explained in de- 
tails in Refiii it is possible to carry out this imaginary- 
time propagation entirely without time discretization er- 
rors for lattice systems. 

To explain how the continuous-time propagation 
works in CTRMC consider first the possible actions of the 
evolution operator on a particular configuration during 
an infinitesimal time step dr. In CTRMC there are two 
possible actions: 1) The configuration changes or, 2) the 
configuration stays unchanged. In order to account for 
the in general non-Markovian nature of the imaginary- 
time evolution operator a weight is also associated to the 
evolving configuration. With this weight it is possible 
to enforce probability conservation for the two actions 
1) and 2) provided the weight is altered if action 2) is 
chosen. This strategy is the same as utilized in so called 
pure (no branching) Diffusion Monte Carlo methods^^*^. 
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The probability of action 1) is given by offdiagonal ele- 
ments of the Hamiltonian and is of the order dr. To void 
the sign problem these elements must all be negative or 
zero. The probability for the configuration to stay un- 
changed, action 2), is p(2) = 1 — p(l) = 1 + Hi C dr, 
where c labels the current configurational state and i are 
other states connected to the current state by an offdiago- 
nal term in the Hamiltonian. This choice of probabilities 
and the form of the infinitesimal time evolution operator, 
1 — H cc dT, implies that the multiplicative weight change 
associated with action 2) is 1 — E L (c)dr where E L (c) 
is the so called local energy, El(c) = H cc + Yli^ic: 01 
configuration c. 

Because p(2) is of the order unity one can simulate the 
continuous-time evolution in the same way as done for 
continuous-time simulations of radioactive decay. One 
generates stochastically a future decay-time Td CC ay ac- 
cording to the distribution e -~{ 1 -p( 2 )) T A^y/ dT and moves 
the configuration directly to this time while multiplying 
its weight with e^ EL ^ Tdocay . At the decay time the type 
of decay, i.e. which state to transfer to, is determined 
dependent on the relative values of the many offdiagonal 
matrix elements leading away from configuration c. In 
this way the simulation is carried out without time-step 
errors. 

An evolving configuration will contribute to observ- 
ables a term proportional to its accumulated weight. 
However it is well known that the weights will be widely 
spread in magnitude^. The sum will thus be dominated 
by just a few terms with the biggest weights thus giving 
essentially only the contributions for a few configurations 
yielding bad statistics. The advantage of the Reptation 
Monte Carlo technique is that this summation over con- 
figuration weights is also carried out stochastically using 
importance sampling based on the weight magnitudes. 
This sampling ensures a much more efficient summa- 
tion as the observables gets relatively many contributions 
from the configurations with the highest weights. 

To make a practical implementation of CTRMC a 
starting configuration is stored as the first element in 
an array of configurations and then propagated a finite 
time interval A T . The resulting configuration is stored 
as the next element in the configuration array while the 
weight associated to the A T propagation is stored as the 
first element of a separate array of weights. This process 
is repeated until the configuration has evolved for a total 
time Ttot = fiA T . The configuration array with n + 1 
elements and the weight array with n elements consti- 
tute then a description of the continuous-time evolution. 
The total weight for this r to t evolution is the product of 
weights in the weight array. The collection of the con- 
figuration and weight arrays will be referred to as the 
polymer. A schematic view of two polymers is shown in 
Fig. [SJ Having filled the polymer arrays the reptation 
move starts. First a random choice of propagation di- 
rection is made: Either the reptation propagation starts 
from the first element of the polymer and propagates 
backwards in time or from the last element propagating 




FIG. 5: Schematic drawing of two polymers or equivalently 
two propagations in state space. Each perpendicular line 
marks an element in the configuration array. Each segment 
between two consecutive perpendicular lines indicates a A T 
time propagation and has an associated entry in the weight 
array. The two polymers illustrate the evolution before(left) 
and after(right) the reptation move. The total weight of these 
two polymers differ only by their first and last segments. The 
dashed segment on the right polymer indicates a removed part 
so that both polymers have the same length. 



forward in time. Depending on the propagation direction 
the appropriate configuration (first or last) is copied and 
propagated for a time interval A T . After this step one es- 
sentially has information about two different evolutions 
that are almost identical except for their first and last 
elements. This is illustrated in Fig. 03 

One now performs a Metropolis accept/reject deci- 
sion based on the relative weights of the two polymers. 
That is: accept the reptation move with probability 
p = min(W /W, 1) where W(W) is the total weight for 
the polymer before(after) the reptation step. This ra- 
tio of weights is easily calculated as it involves only the 
weights at the ends of the polymers. It depends clearly 
on the value of A T . We have found it most efficient to 
adjust A r such that the resulting acceptance probability 
is about 1/2. 

Diagonal observables can easily be extracted from any 
stored configuration in the polymer. By extracting ob- 
servables from the middle block in the polymer we en- 
sure that observables are picked according to the forward- 
walking procedure2i with a forward propagation time of 

Ttot/2. 

The only adjustable parameter in this scheme is r to t 
which ideally should be as long as possible. In this article 
we have taken it to be ~ 40 J -1 , and 80 J -1 close to the 
RFC-point. 

CTRMC can, as any other projection Monte Carlo 
techniques, be improved substantially by using a guid- 
ing function^. The guiding function should be as close 
as possible to the true ground state wave function. In 
our simulations we have chosen a guiding wave function 
that is biased towards having many plaquettes with par- 
allel dimers. We optimize the guiding wave function by 
performing a small trial run before the actual run. In 
the trial run the guiding wave function is optimized in 
such a way as to yield a minimal variance of the lo- 
cal energjiS. At the RK-point, where the ground state 
is explicitly known, the resulting quantum Monte Carlo 
simulation using the exact ground state as guiding wave 
function reduces to a classical Monte Carlo simulation 
which also can be employed to measure properties of ex- 
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ited statea 17 ! 24 . 25 ! 26 . 

There is another class of Projector Monte Carlo meth- 
ods where the non-Markovian character of the evolution 
operator is taken care of by introducing extra processes 
such as replication and decimation of copies of the sys- 
tem sometimes referred to as walkers^Z* 2 ^. Having such 
a changing number of walkers is undesirable as it will 
need some form of population control to prevent fluctu- 
ations from killing all walkers or filling the entire com- 
puter memory. Population control leads again to a sys- 
tematic bias of the results which must be corrected for 
by reweighingSii. In CTRMC there is only one walker, 
thus these problems are avoided. It is also possible to get 
around these problems using a population with a stochas- 
tic reconfiguration of a constant population of walkers as 
demonstrated in Refi^S. 



IV. PLAQUETTE PHASE 
A. Restoration of rotational symmetry 

The columnar ordering phase is in part distinguished 
from the plaquette phase by having a preferred dimer ori- 
entation. Thus a suitable order parameter distinguishing 
these phases is one that detects any orientational pref- 
erences of the dimers. Such an order parameter can be 
constructed by finding an operator that changes sign un- 
der a 7r/2 lattice rotation and a subsequent translation. 
The translation is needed to make the plaquette state 
invariant under this transformation. 

Leung et al& proposed to detect this by measuring the 
difference between the number of vertical and horizontal 
dimers 

M vh = -^(N v - N h ) . (3) 

Here N v (Nhh) is the total number of vertical (horizontal) 
dimers. In the ideal columnar state M 2 h = 1 whereas in 
the ideal plaquette state M 2 h = A/L 2 which vanishes in 
the thermodynamic limit. M v h is therefore suitable as 
an order parameter. 

Leung et al. studied M v h and its Binder ratio up to 
system sizes L — 8 and concluded that there is a phase 
transition at V c = —0.2. However the crossing points of 
the Binder ratios vary significantly with system size so 
bigger systems are needed to determine the location of 
the phase transition more accurately. Progress on bigger 
system sizes, up to L = 20 were obtained in Refitt using a 
Diffusion Monte Carlo method, but no finite-size analysis 
of the data leading to a definite estimate for the location 
of the phase transition were given. 

In order to find the location of the phase transition we 
have extracted (M 2 ^) 1 ' 2 from our simulations for differ- 
ent system sizes up to L = 32. In Fig. [5] we have plotted 
the results for \J (M v h) 2 as a function of L for different 
values of V. The solid lines are best fits to the func- 
tional form a/L 2 . In the inset we have plotted 




FIG. 6: Color online. The order parameter (M^ h } 1/2 vs. l/L 
for different values of V. The top curve has V = and in- 
creases in steps of 0.1 for lower curves. The solid lines are fits 
to the form y>H- a/L 2 . The inset shows the values of fi vs. 
V. The inset contains data for more values of V than shown 
in the main figure. 



the resulting infinite size extrapolation value (/i) where 
also results from other values of V is included. From 
this figure we get our best estimate for the location of 
the phase transition between the columnar and plaque- 
tte state: V c = 0.60 ± 0.05. 

We have also extracted higher moments of the M v h 
distribution. In Fig. [7| we show the Binder ratio 
(-^Ch) / (Mvh) 2 as functions of V for different system sizes. 
For large values of V all curves have values close to 3, 
consistent with a Gaussian distribution with zero mean. 
The curves for different system sizes do not cross in a 
single point. However, as seen in the inset of Fig. [7] the 
crossing points for curves of system sizes L and L + 4 
seems to converge to a point close to V c = 0.6 consistent 
with the estimate above. However the rather large error 
bars on the crossing points do not constrain this estimate 
further. 

The phase transition between the columnar and pla- 
quette phases on the square lattice is often conjectured 
to be a first order transition. This is consistent with what 
happens on the hexagonal lattice^ and in other models 
with similar phase transitions^*^. The limited system 
sizes studied here makes it however rather difficult to 
verify this conjecture. We do not observe any hints of 
discontinuities in the differentiated ground state energy 
for the largest system sizes. Neither do we see any dis- 
continuities in the way the order parameter approaches 
zero. However this does not rule out a weak first order 
transition that might only be visible in simulations using 
larger lattices. 
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FIG. 7: Color online. Binder ratios vs. V for different system 
sizes. The top curve on the left side corresponds to linear 
system size L = 8. L increases successively by 4 for lower 
curves (on the left site). The inset marks the crossing points 
for L and L + 4 as functions of L. 

B. Breaking of translational symmetry 

The ideal columnar state in Fig. ^ is invariant under 
translation by one lattice spacing perpendicular to the 
orientation direction of the dimers. The plaquette phase 
breaks this symmetry. One could imagine an intermedi- 
ate phase that breaks both this translational symmetry 
as well as the rotational symmetry discussed in the previ- 
ous section. A phase like this could resemble that shown 
in the upper left panel of Fig. 0] and be described for in- 
stance as the product state in Eq. [21 but with different 
amplitudes for the vertical and horizontal directions. 

An order parameter which detects the breaking of 
translational symmetry by one lattice spacing in the di- 
rection perpendicular to the orientation of the majority 
of dimers is 

M trans = A^[0( Mvh )iVvv(r)(-ir 

plaq 

+ 6(~M vh )N hh (r)(~iyy] (4) 

where iVhh {N vv ) takes the value one for a plaquette 
with two horizontal (vertical) dimers and zero otherwise. 
#(M v h) = 1 if there are more vertical than horizontal 
dimers and otherwise, 0(—x) = 1 — 6(x). r denotes the 
midpoint coordinate of a plaquette; its component are 
integers in units of the lattice spacing. M trans is invari- 
ant under n/2 lattice rotations and changes sign when a 
configuration is translated one lattice spacing along the 
direction perpendicular to the orientation of the majority 
of dimers. Thus A/trans — in the columnar phase, while 
it takes the value ±1/2 for the ideal plaquette states. The 




FIG. 8: Color online. Main panel shows the transla- 
tional symmetry breaking order parameter vs. inverse lin- 
ear system size for different values of V. The values are 
for the curves from top to bottom on the right side V = 
0.8, 0.6, 0.4, 0.2, 0.0, -0.2, -0.4, -1.0, -2.0. The inset shows 
the same order parameter, but plotted as a function of V for 
different system sizes. Here the smallest system size is for the 
top-most curve. 



inset of Fig. |H1 shows (M t 2 rans ) vs. V for different system 
sizes. For a fixed system size this order parameter stays 
small for low V and increases towards a maximum at 
about V = 0.6 and then decreases again towards V = 1. 
The decrease at high values of V is partly due to the 
reduced number of flippable plaquettes. For all the sys- 
tem sizes investigated the order parameter decreases with 
system size. 

In order to see the finite-size behavior of this decrease 
we have plotted the order parameter as a function of in- 
verse linear system size in the main panel of Fig. |H1 It 
appears that the order parameter extrapolates to for 

V < 0.2. Based on this we rule out translational sym- 
metry breaking for V < 0.2. A close inspection of the 

V = 0.4 curve reveals a downward curvature similar to 
that seen more clearly for the V = 0.2 curve, thus we 
believe that the V = 0.4 will also extrapolate to zero. 
In contrast the V = 0.6 curve reveals an upward cur- 
vature indicating an extrapolation to a finite value and 
therefore breaking translational symmetry. The V = 0.8 
curve lies lower than the V = 0.6 but has also a slight 
upward curvature. Thus we conclude that the transla- 
tional symmetry breaking happens for V ~ 0.4 — 0.6. 
However larger system sizes are necessary to determine 
a more precise value. 

A critical value of V ~ 0.4 — 0.6 is close to the value 
of V where rotational symmetry gets restored as found 
in the previous section. Thus if a phase that breaks both 
translational and rotational symmetry exists, it can only 
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FIG. 9: Average value of flip operator for different values of 
V. 

do so in a rather narrow region of V close to the rotational 
symmetry breaking transition. 
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FIG. 10: Offdiagonal correlation function (-Fb-Fi) as function 
of diagonal separation r; = at V = 0.9 on a L = 32 

lattice. The inset shows the behavior of for different 

values of L. The line in the inset is the best linear fit to 
the largest system sizes (excluding L = 8) and is A^/2 = 
0.023 + 0.17/L. 



V. OFFDIAGONAL CORRELATION 
FUNCTIONS 

The resonating plaquette phase is most directly probed 
by looking for resonating dimers on one of the four sub- 
lattices. A resonating plaquette is characterized by a 
finite expectation value for the off-diagonal flip operator 

F i = 1^(111 + 111)^1, (5) 

where i labels the plaquette. To get an impression of 
how much the dimers are resonating we plot in Fig. Othe 
average value of the flip operator per plaquette. This can 
be calculated directly from measuring the ground state 
energy per plaquette E and the potential energy 

±- 2 Y,{Fi) = -E + V{^), (6) 

i 

Nf is the number of plaquettes with parallel dimers in a 
given configuration. We see that the dimers resonate the 
most at V = 0. 

To get a picture of how the resonating dimer plaquettes 
are correlated on the lattice we measure the correlation 
function (FiFj). In the plaquette phase this correlation 
function should show long range order when i and j both 
belong to the sublattice with resonating dimers. In order 
to measure {FiFj) we used the Feynman-Hellman theo- 
rem. An operator —aFiFj was added to the Hamiltonian 
and the ground state energy was measured as a function 
of a. Taking small values of a the expectation value of 
FiFj was determined as (FiFj) = — g^~\a^o- 



In Fig. UUIwe have plotted (F Fi) for V = 0.9 at di- 
agonal spatial separations: i denotes a plaquette with 
coordinates (i,i). The zero result for i — 1 is an ex- 
act result as two diagonally adjacent plaquettes cannot 
both be flippable. The curve exhibits clear oscillations in 
the bulk indicative of a resonating plaquette phase. To 
find out if these oscillations are also present in the ther- 
modynamic limit, we have defined a quantity A^/ 2 that 
measures the magnitude of the oscillations in the bulk: 
A L/2 = (F F L/2 ) - {FoFl/^). In the inset of Fig. QUI 
we have plotted function of 1/L. A L / 2 ex- 

trapolates to a finite value in the thermodynamic limit, 
consistent with the presence of a resonating plaquette 
phase. One should however note that the magnitude of 
these bulk oscillations are small compared to what is ex- 
pected in the ideal plaquette state, Eq.|2 A Aoo = 0.023 
for V — 0.9 corresponds to 9% of the value expected for 
one of the ideal plaquette states. 

To contrast this finding at V = 0.9 we have repeated 
the same measurements at the RK-point. The obtained 
result for (FoFi) at the RK-point is shown in Fig. ^] 
for an L = 48 lattice. Although visible the bulk oscilla- 
tions are much smaller in this case. In fact they vanish 
completely in the thermodynamic limit as can be seen 
from the inset of Fig. llll which shows the finite-size scal- 
ing of A L / 2 which vanishes as a power law L~ 9 . The 
bulk decay at the RK-point can be calculated analyt- 
ically by calculating the correlation function for having 
two parallel dimers on a plaquette displaced from another 
plaquette also with two parallel dimers. Using the Pfaf- 
fian technique^ and the Green function given in Ref^ it 
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FIG. 11: Offdiagonal correlation function (FoFi) as function 
of diagonal separation fi = (i, i) at the RK-point on an L — 48 
lattice. The inset shows A L / 2 vs - l/L. The line is the best fit 
to a power law A L/2 ~ l/L 3 with g = 2.02 ± 0.03. 
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FIG. 12: Color online. <M c 2 ol ) vs. l/L for various val- 
ues of V. From top to bottom the values of V are V = 
0.960, 0.970, 0.975, 0.980, 0.985, 0.990, 0.995, 0.998, 1.000. The 
dashed lines are guides to the eye. 



is easy to show that asymptotically (FoFi) ~ (—iy/i 2 ^ 
thus the exact value of g = 2. Our numerical result is 
consistent with this. 



VI. TRANSITION TO RVB LIQUID 



The columnar order parameter 

$>h(r)(-ir 



Ml 



1 

iZ 1 



, plaq 



+ ($>v(r)(-iP 

■ plaq 



(7) 



where N^r) (N v (f)) is the number of horizontal (ver- 
tical) dimers surrounding the plaquette at r, was pro- 
posed in RefpiS as a mean to detect the columnar order, 
l-^coil = 1/2 in the ideal columnar state. The columnar 
order parameter is constructed so that it is for any state 
invariant under lattice rotations. The plaquette state is 
not invariant under lattice rotations alone. It is invariant 
under the combined operation of a rotation and a trans- 
lation, thus in fact the columnar order parameter for the 
ideal plaquette state is |M co i| = 1/V8 in the thermo- 
dynamic limit. Therefore the columnar order parameter 
does not distinguish between the columnar and plaque- 
tte phases in a useful way. However the columnar order 
parameter is zero in the RVB liquid and in the staggered 
state, thus the phase transition from the plaquette phase 



to the RVB liquid can be detected using the columnar 
order parameter. 

In Fig.^|we show the size dependence of M^ ol for dif- 
ferent values of V close to the RK-point. For all curves 
the order parameter decreases with increasing system 
size. As V is moved away from 1 the decrease becomes 
slower and appears to saturate to a finite value, at least 
for the curves with V < 0.99. This saturation to a finite 
value is not manifest for the curves with V close to 1 for 
the system sizes considered here. 

In order to search for a possible phase transition at 
V < 1 we show the Binder ratio (M^ ol ) / (M^ ol ) 2 as func- 
tions of V for different system sizes in Fig.^H The curves 
for different sizes do not cross in a single point. There is 
a tendency that the crossing points of curves for nearby 
system sizes move towards 1 as the system size increases. 
Interpreted this way we conclude that there is no evi- 
dence for a phase transition to the RVB liquid at V < 1 
and that there is significant finite-size corrections to scal- 
ing even for the largest system sizes considered here. 

To check the columnar order parameter results at the 
RK-point we have in addition employed a very effective 
directed-loop Monte Carlo method^ which only is appli- 
cable exactly at the RK-point. Fig. ^] shows the re- 
sults plotted so as to expose the logarithm present in 
the leading asymptotically result (M^ ol ) = C\og(L)/L 2 . 
This result can be calculated analytically using the Pfaf- 
fian technique with the Green function given in Ref. 33 . 
The directed-loop Monte Carlo results agree with this 
asymptotic behavior, see Fig. ^] both when all topolog- 
ical sectors are included in the sampling (lower curve) 
and when the sampling is restricted to the zero winding 
number sector (upper curve). We find that the value of 
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V 

FIG. 13: Color online. The Binder ratio (A4 4 ol )/{M c 2 ol ) 2 vs. 
V for different system sizes. From top to bottom on the left 
side the curves are for sizes L = 16, 24, 32 and 48. 

C = 0.63 ± 0.01. Fig. [H also shows the CTRMC data 
from Fig. El (squares). They coincide with the directed- 
loop Monte Carlo data restricted to the zero winding 
number topological sector. 

VII. CONCLUSION 

We have applied a continuous-time variant of the rep- 
tation quantum Monte Carlo method to study the square 
lattice QDM. In particular we found the location of the 
phase transition between the columnar and plaquette 
phase to be at V c = 0.60 ± 0.05. This phase transition 
happens at a positive V c and excludes the resonating pla- 
quette state as the ground state when the Hamiltonian 
consists of just the kinetic term. 

The estimate for the location of the phase transition is 
based on the finite-size extrapolation of the order param- 
eter and the apparent convergence of the Binder ratios. 
Thus there is a possibility that our estimate is strictly a 
lower bound as we cannot rule out a scenario were the or- 
der parameter is very small but finite, and where a phase 
transition is only seen for very large lattices. Our results 
shine little light on the thermodynamic properties of the 
phase transition. The crossing points for the Binder ra- 
tio of the rotational symmetry breaking order parameter 
converge rather slowly, thus there are significant finite- 
size corrections to scaling if the transition is continuous. 
If the transition is first order, it is rather weakly so, as we 
see no evidences for discontinuities in the differentiated 
ground state energy nor in the order parameter for the 
system sizes considered here. 

That finite-size effects are very important in these 
studies is also apparent from the average dimer densi- 
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FIG. 14: Directed-loop Monte Carlo simulation results show- 
ing (M 2 ol )L 2 vs. L at the RK-point on a semi-log plot. In 
the upper curve the data points are obtained by restricting 
the sampling to the zero winding number sector. The squares 
are centered about the Reptation Monte Carlo results gotten 
from Fig. 1121 In the lower curve the data points are col- 
lected from all topological sectors. The lines are best fits to 
the data using the functional form (Af 2 ol ) — C\og(L)/L 2 and 
gives within error bars the same C = 0.63 ± 0.01 for both 
curves. The largest system size simulated had L = 2048. 

ties in open boundary geometries shown in Fig. 0] where 
one can clearly see plaquette patterns both far into the 
columnar phase and at the RK-point. 

We have also considered the possibility of an intermedi- 
ate phase in between the columnar and plaquette phase 
breaking both rotational and translational symmetries. 
In order to do so we have constructed an order parameter 
that is insensitive to rotations, but detects the breaking 
of translational symmetry in the direction perpendicular 
to the orientation of the majority of dimers. We find a 
phase transition occurring at roughly the same value of 
V as the rotational symmetry gets restored. Thus if such 
an intermediate phase exists it is confined to a narrow 
region in phase space close to V ~ 0.6. 

To strengthen our argument for the existence of the 
plaquette phase we have also measured the off-diagonal 
dimer-fiip correlation function inside the plaquette phase 
and performed a finite-size scaling of the results. This re- 
veals long-range (staggered) order in the thermodynamic 
limit consistent with the existence of a plaquette phase 
with resonating dimers. However the strength of the pla- 
quette pattern seen at V = 0.9 is rather weak. It is 
only about 9% of the ideal plaquette state value. At the 
RK-point we find that the plaquette pattern vanishes as 
the square of the linear system size in agreement with 
analytic calculations. 

We have also considered the transition from the pla- 
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quette phase to the RVB liquid as measured by the 
columnar order parameter. To verify the CTRMC data 
for the columnar order parameter at the RK-point we 
measured the columnar order parameter using a directed- 
loop Monte Carlo algorithm. These methods gave the 
same results, and the directed-loop Monte Carlo algo- 
rithm was employed to verify the analytic prediction for 
the finite-size scaling of the columnar order parameter up 
to linear system sizes L = 2048. 

From the columnar order parameter data away from 
the RK-point, there are no evidences for a phase tran- 
sition occurring for V < 1. The crossing points of the 
Binder ratios for different system sizes from L = 16 to 
L = 48 move towards higher values of V as the sys- 
tem sizes are increased. Thus we conclude that there is 
significant finite-size corrections even at the largest sys- 
tem sizes considered here. Large finite-size effects is to 
be expected as the effective height-model describing the 



system close to the RK-point contains a dangerous irrel- 
evant operatorJ^i, making the extraction of scaling pa- 
rameters from finite sized samples a complicated issued. 
More numerical work on larger lattices combined with a 
proper finite-size scaling ansatz is needed to extract the 
proper critical behavior of this phase transition. 
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